The simulation problem below presents a vertical inverted pendulum with an excited base. This problem becomes interesting due to the fact that given certain initial angles from vertical and the right amplitudes and frequency's of excitation, the pendulum stays upright without falling over. Simulations provided below model this behavior and display the proper initial conditions that promote stabilization of the inverted mass. Animations are created to display the inverted pendulum stabilizing across variations in initial conditions.
Below can be seen two images that setup the problem we are looking at. The diagram of the model can be seen with the excited base moving up the Y axis and inverted mass connected with a rod free to move in the X-Y plane. The mass is only influenced by the force of gravity and restricted by the connection to the exicted base. The angle $\theta$ is defined as positive from the vertical Y-axis moving clockwise in the diagram.

The bond graph below shows the physical model of the system and how each part is interacting such that we can derive equations to solve for the motion of the inverted mass.

Velocity equations for X, Y, and $\theta$ directions
$$v_y = \dot{y} - l \dot{\theta} \sin(\theta) $$$$v_x = l \dot{\theta} \cos(\theta) $$$$ \dot{\theta} = \frac{p_x}{l\cos(\theta)\ m} $$Momentum equations
$$p_x = m l \dot{\theta} \cos(\theta) $$$$p_y = m\left[\dot{y} - \frac{\sin(\theta)\ p_x}{\cos(\theta)\, m}\right] $$Derivatives of momentum equations $$\dot{p_y} = m\left[\ddot{y} - \frac{\sin(\theta)\ \dot{p_x}}{\cos(\theta)\ m} - \frac{p_x\ \dot{\theta}}{m\ \cos^2{\theta}} \right] $$ $$\dot{p_x} = m g \sin(\theta)\cos(\theta) + m \ddot{y}\sin(\theta)\cos(\theta) - \frac{\sin(\theta)}{\cos(\theta)} \dot{\theta} p_x$$
Since we plan to do several simulations and visualize them, in this section we define several functions that allow us to simplify how this is done.
To solve this non-linear system of ordinary differential equations, we need to start from an initial state, numerically integrate along the state derivatives to find the next state, and repeat as desired
For integration, we will use scipy.integrate's solve_ivp(), which takes in a range of time (t_span) through which to integrate, a function handle (func) that takes in the time and state and returns the derivatives, and the initial condition (initial) of the state.
Since we wish to solve this problem using several different parameters, we will define a function get_func() that returns a function handle that we can pass to solve_ivp(). We also return a function that includes additional state information.
from copy import deepcopy
import numpy as np
def get_func(params: dict[str, float]):
# Make a copy of the input parameters to ensure the returned function
# does not change when the original parameters dictionary changes
params = deepcopy(params)
# Extract parameters from dict
theta = params["theta"]
ampl = params["amplitude"]
freq = params["frequency"]*2*np.pi # Convert from Hz to rad/s
leng = params["l_arm"]
mass = params["m_mass"]
g = params["g"]
# Get initial condition (add position into the system so solve_ivp can integrate it for us)
px = 0
initial = [px, theta]
# Create derivative calculator function to be used by solve_ivp()
# state = [px, theta]
def func(t, state):
state = deepcopy(state)
px = state[0]
theta = state[1]
# Sine and cosine of theta, for computational speed
sin_theta = np.sin(theta)
cos_theta = np.cos(theta)
tan_theta = np.tan(theta)
# Position of cart and its derivatives
X = 0
Y = ampl*np.sin(freq*t)
d_Y = ampl*freq*np.cos(freq*t)
dd_Y = -(freq**2)*Y
# State derivatives
d_theta = px/(leng*cos_theta*mass)
d_px = mass*g*sin_theta*cos_theta + mass*dd_Y*sin_theta*cos_theta - tan_theta*d_theta*px
# Other state variables
x = X + leng*np.sin(theta) # x-position of pendulum mass
y = Y + leng*np.cos(theta) # y-position of pendulum mass
py = mass*(d_Y - tan_theta*px/mass) # y-momentum of pendulum mass
d_py = mass*(dd_Y - tan_theta*d_px/mass - px*d_theta/(mass*cos_theta**2)) # first time-derivative of y-momentum of pendulum mass
# Concatenate state derivatives
d_state = [d_px, d_theta]
# Create current state dict
state = {}
state["t"] = t
state["X"] = X
state["Y"] = Y
state["x"] = x
state["y"] = y
state["px"] = px
state["py"] = py
state["theta"] = theta
state["d_Y"] = d_Y
state["d_px"] = d_px
state["d_py"] = d_py
state["d_theta"] = d_theta
state["dd_Y"] = dd_Y
state["tan_theta"] = tan_theta
return d_state, state
def func_wrap(t, state):
'''
Since solve_ivp() needs a function that only returns the state
derivatives, we create a wrapper for func() that discards the rest of
the state and returns only the state derivatives
'''
d_state, state = func(t, state)
return d_state
return func, func_wrap, initial
solve_ivp()¶Now that we have a function that calculates the state derivatives given the state, we need to propagate the solution throughout time. The function solve_problem() was designed to use the get_func() function we defined earlier, calculate the solution using scipy.integrate's solve_ivp(), and consolidate the interesting parts of the solution and the parameters given to define it.
from scipy.integrate import solve_ivp
import pandas as pd
import numpy as np
def solve_problem(
params: dict[str, float],
t_eval: np.ndarray,
name: str,
rtol: float = 1e-12,
atol: float = 1e-12,
):
func, func_wrap, initial = get_func(params)
output = solve_ivp(
func_wrap,
(min(t_eval), max(t_eval)),
initial,
t_eval=t_eval,
method="RK45",
rtol=rtol,
atol=atol,
)
ts = output.t
ys = output.y
states: list[dict[str, float]] = []
for t, y in zip(ts, ys.T):
d_state, state = func(t, y)
states.append(state)
df = pd.DataFrame(states)
solution = {
"params": params.copy(),
"data": df,
"name": name,
}
return solution
To visualize the solutions, we can create a static plot or an animation of the pendulum. Since we have several solutions we would like to plot together, we create functions to consolidate the code. plot_time() and plot_ani() each take in a list of dictionaries outputted as the solution from solve_problem().
from typing import Any
from matplotlib.lines import Line2D
import matplotlib.pyplot as plt
from matplotlib import animation
import numpy as np
import pandas as pd
def plot_time(solutions):
fig = plt.figure()
ax = fig.add_subplot()
for solution in solutions:
name = solution["name"]
df = solution["data"]
t_vals = df.get("t").to_numpy()
theta_vals = df.get("theta").to_numpy() * 180 / np.pi
ax.plot(t_vals, theta_vals, label=f"theta: {name}")
plt.legend()
plt.title("Theta vs. Time")
plt.xlabel("time (s)")
plt.ylabel("theta (deg)")
return fig
def plot_ani(solutions: list[dict[str, Any]], interval: int = 10):
fig = plt.figure()
ax = fig.add_subplot()
all_x = np.array([])
all_y = np.array([])
# Initialize arrays to store the animation plots and data to use to update them
num_solutions = len(solutions)
plots: list[dict[str, Line2D]] = [{} for n in range(num_solutions)]
data: list[dict[str, np.ndarray]] = [{} for n in range(num_solutions)]
for i, solution in enumerate(solutions):
name: str = solution["name"]
df: pd.DataFrame = solution["data"]
data[i]["t_vals"] = df.get("t").to_numpy()
data[i]["X_vals"] = df.get("X").to_numpy()
data[i]["Y_vals"] = df.get("Y").to_numpy()
data[i]["x_vals"] = df.get("x").to_numpy()
data[i]["y_vals"] = df.get("y").to_numpy()
(plots[i]["cart"],) = ax.plot([], [], "o", c="b", label=f"cart: {name}")
(plots[i]["line"],) = ax.plot([], [], "-", c="k")
(plots[i]["pend"],) = ax.plot([], [], "o", label=f"pendulum: {name}")
num_frames = len(data[i]["t_vals"])
all_x = np.concatenate((all_x, data[i]["X_vals"].copy(), data[i]["x_vals"].copy()))
all_y = np.concatenate((all_y, data[i]["Y_vals"].copy(), data[i]["y_vals"].copy()))
def update_points(n):
returns = []
for i in range(num_solutions):
plots[i]["cart"].set_data(([data[i]["X_vals"][n]], [data[i]["Y_vals"][n]]))
plots[i]["line"].set_data(
(
[data[i]["X_vals"][n], data[i]["x_vals"][n]],
[data[i]["Y_vals"][n], data[i]["y_vals"][n]],
)
)
plots[i]["pend"].set_data(([data[i]["x_vals"][n]], [data[i]["y_vals"][n]]))
returns.extend([plots[i]["cart"], plots[i]["line"], plots[i]["pend"]])
returns = tuple(returns)
return (*returns,)
ani = animation.FuncAnimation(
fig, update_points, num_frames, interval=interval, blit=True, repeat=True
)
min_x, max_x = min(all_x), max(all_x)
min_y, max_y = min(all_y), max(all_y)
plt.xlim(min_x - 0.5, max_x + 0.5)
plt.ylim(min_y - 0.5, max_y + 0.5)
plt.legend()
plt.title("Pendulum Animation")
plt.xlabel("x (m)")
plt.ylabel("y (m)")
return ani
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams["animation.html"] = "jshtml"
# -----------------------------
# Global Parameters
# -----------------------------
params: dict[str, float] = {}
params["l_arm"] = 1.0 # length of pendulum arm [m]
params["m_mass"] = 1.0 # mass of pendulum mass [kg]
params["g"] = 9.8 # gravitational force [m/s^2]
# Time
t_start: float = 0
t_end: float = 25
t_increment: float = 0.05 # visualization time step duration (sec)
# NOTE: playback speed is t_increment/(interval/1000)
t_span = (t_start, t_end)
t_eval = np.arange(min(t_span), max(t_span)+t_increment, t_increment)
params["amplitude"] = 0
params["frequency"] = 0
params["theta"] = np.deg2rad(60)
solution = solve_problem(params, t_eval, "validation")
fig = plot_time([solution])
As we can see from the validation case, the pendulum is behaving as expected. It does not appear to have any unknown forces acting on it, and it has a minimum and maximum angle of +/- 60°, which is where it was released.
Here we simulate the inverted pendulum with the oscillating cart at the recommended amplitude and frequency.
# Time
t_start: float = 0
t_end: float = 5
t_increment: float = 0.005 # visualization time step duration (sec)
# NOTE: playback speed is t_increment/(interval/1000)
t_span = (t_start, t_end)
t_eval = np.arange(min(t_span), max(t_span)+t_increment, t_increment)
# Run simulations
solutions = []
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(10)
solutions.append(solve_problem(params, t_eval, "10°"))
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(20)
solutions.append(solve_problem(params, t_eval, "20°"))
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(30)
solutions.append(solve_problem(params, t_eval, "30°"))
fig = plot_time(solutions)
We see that the for the three initial angles, 10°, 20°, and 30° away from vertical, the inverted pendulum is stable with amplitude 0.1 meter and frequency 20 Hz.
Although not necessary, it is interesting to see where the 0.1 meter, 20 Hz oscillation stops working. Adding in a 40° initial angle to the simulation shows that the oscillation no longer works at this initial angle. Here, we see the 40° pendulum begins to rotate continuously around the cart, as the cart maintains a higher level of energy in the system than it starts with.
In the case of 10°, 20°, and 30°, we see a sort of "destructive interference" between the oscillation of the pendulum and the oscillation of the cart. However, when we get to 40°, we see "constructive interference". This implies that there would be a better frequency/amplitude for a pendulum starting at this initial angle.
# Time
t_start: float = 0
t_end: float = 1.5
t_increment: float = 0.005 # visualization time step duration (sec)
# NOTE: playback speed is t_increment/(interval/1000)
t_span = (t_start, t_end)
t_eval = np.arange(min(t_span), max(t_span)+t_increment, t_increment)
# Run simulations
solutions = []
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(10)
solutions.append(solve_problem(params, t_eval, "10°"))
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(20)
solutions.append(solve_problem(params, t_eval, "20°"))
params["amplitude"] = 0.1
params["frequency"] = 20
params["theta"] = np.deg2rad(30)
solutions.append(solve_problem(params, t_eval, "30°"))
params["amplitude"] = 50
params["frequency"] = 0.1
params["theta"] = np.deg2rad(40)
solutions.append(solve_problem(params, t_eval, "40°"))
fig = plot_time(solutions)
Just for fun, we can use our animate function to show a time-based animation of all four pendulums. Before starting it, note the starting positions are prescribed, and then see the 40° pendulum make full rotations about the cart.
from IPython.display import HTML
# Plot solutions
ani = plot_ani(solutions, interval = 15)
plt.axis('equal')
HTML(ani.to_jshtml())